---
title: "Replication of Ihme & Tausendpfund (2018) - Stage 1 and 2"
author: "Flávio Azevedo, Deliah Bolesta, Leticia Micheli"
date: "23/10/2022"
output:
  html_document:
    toc: true
    toc_float: true
  pdf_document: default
editor_options: 
  chunk_output_type: console
---

```{r libraries and functions, echo=F, message=FALSE, warning=FALSE}
## Libraries needed for the script to run ----
#install.packages("see", repos = "https://easystats.r-universe.dev")
library(apaTables) # creates APA tables to Word
library(BayesFactor)
library(car) # for ANOVA
library(curl)
library(difR)
library(effects)
library(emmeans) # for estimated marginal means
library(extrafont)
font_import()
loadfonts(device = "win")
library(ggplot2) # for data visualization
library(ggpubr) # for adding a stats layer
library(grid)
library(gridExtra)
library(ggstatsplot)
library(ggthemes)
library(gmodels)
library(haven) # for codebook
library(knitr)
library(labelled) # for codebook
library(latticeExtra)
library(ltm) # for Cronbach's Alpha
library(MASS)
library(mirt)
#library(pairwiseComparisons) # pairwise comparisons
library(pwr)
library(sjstats) # for effect sizes in ANOVA
library(stringr) # for the detection of patterns in strings
library(tidyr)
library(wesanderson)
library(dplyr)
options(digits = 3)

## Functions ----
# z Test
z.test <- function(x1, x2, sigma1, sigma2){
  z <- (x1 - x2)/(sqrt((sigma1)^2 + (sigma2)^2))
  print(z)
}

# centering with "scale"
center_scale <- function(x) {
    scale(x, scale = FALSE)
}
```

# Power calculations

This is the power analysis used to calculate the sample size for verification attempts made on Ihme and Tausendpfund (2018). It is based on information extracted directly from the paper. 

```{r record-variables, include = FALSE}
#These variables come directly out of the database. (In the future, a script will retrieve them.) If any are blank/not present, use the default empty string. Make sure the text is surrounded by quotes, or the script won't work. 

pdf_filename = 'Ihme_JournExpPoliSci_2018_xYbO' #Change this to your code!
coded_claim3b = 'The study conducted a 2 × 2 × 3 ANCOVA for the dependent variable (political knowledge test score) using political interest as the covariate. The independent variables were the participants’ gender (male vs. female), the participants’ ﬁeld of study (psychology vs. politics), and stereotype activation (stereotype not activated vs. stereotype activated by gender-question vs. stereotype activated by gender difference statement). Subjects were randomly assigned to the respective gender stereotype conditions.'
coded_claim4 = 'Controlling for political interest, the interaction of stereotype activation and gender (F(2,364) = 6.17, p = 0.002, η² = 0.03) was significant.'
coder = 'Andrew Tyner'
original_coded_effect_size = '0.03'
original_coded_sample_size = '377 students'
```

## Variable standardization

The variables needed for the power analysis are first extracted into the necessary formats. 

```{r convert_variables, include = TRUE}
original_test_statistic_type = 'f'
original_test_statistic_numeric = 6.17
original_effect_size_type = 'eta squared'
original_effect_size_numeric = .03
original_sample_size_numeric = 377
```

Many papers do not report an effect size, and where they are present, there may be ambiguity about which effect size is being reported, or conflictings between the reported test statistic and effect size. We therefore re-calculate all effect sizes, and resolve any disagreements in favor of the reported test statistic (e.g. 'trusting' the papers report of `t` over `d`, if both were reported.)

```{r recalculate_effectsize}
original_recalculated_effect_size_numeric = 0.0327895 #Note that this variable should have the same units as original_effect_size_numeric
```

## Power analysis calculation

For all papers in SCORE, there are two power calculations used to derive the sample sizes that will be used for RR attempts. The first is 90% power, assuming that the true effect size is 75% of that reported in the original study. The second is 90% power, assuming that the true effect size is 50% of that reported in the original study. 

For this paper, these calculations were conducted following the method in `SCORE power scripts.R` for ANCOVA, shown below: 

```{r actual-power-analysis, include = TRUE}
# Remember to use set-seed if this involves any bootstrapping or other processes with randomization!

Fvalue  <- 6.17
df1 <- 2
df2 <- 364
  
p_eta2 <- (Fvalue * df1)/((Fvalue * df1) + df2)
cohenf2 <- p_eta2/(1 - p_eta2)

library(pwr)
SSRP_75 <- ceiling(pwr.f2.test(u = df1, v = NULL, f2 = cohenf2 * .75^2, sig.level = .05, power = .9)$v) + df1 + 1
SSRP_50 <- ceiling(pwr.f2.test(u = df1, v = NULL, f2 = cohenf2 * .5^2, sig.level = .05, power = .9)$v) + df1 + 1

rr_stage1_sample_size = SSRP_75
rr_stage2_sample_size = SSRP_50
```

## Sample size report

For this paper, replication/reproduction will use **a sample of `r rr_stage1_sample_size` for stage 1 data collection**, increasing to **a total of `r rr_stage2_sample_size`** if stage 2 data collection is conducted. 

## Data validation and saving. 

Finally, here is some code that checks that we did indeed get all the necessary values for this power analysis, and then saves those values. This section is the same for all power analyses in SCORE. 

```{r check-data, include = FALSE}
pdf_filename
original_test_statistic_type
original_test_statistic_numeric
original_effect_size_type
original_effect_size_numeric
original_recalculated_effect_size_numeric
original_sample_size_numeric
rr_stage1_sample_size
rr_stage2_sample_size
```

# Data Stage 1

```{r stage 1 read data, echo=F, message=FALSE, warning=FALSE}
# Data from Stage 1
data_stage1 <- sjlabelled::read_spss("analytical_data_stage1.sav")
# Original data from Ihme & Tausendpfund
originaldata <- haven::read_spss("original_data.sav")
# authors included a filtr variable (Filter) that assigned a value of zero to anyone who
# were neither female nor male, did not intend to answer honestly and who weren't German
originaldata <- originaldata[which(originaldata$Filter == "1"),]
```

```{r stage 1 exclude cases, echo=F, message=FALSE, warning=FALSE}
## data_stage1 exclusion ----
data_stage1 <- data_stage1[which(data_stage1$Complete == "Yes"),]
```

```{r stage 1 data_stage1 wrangling, echo=F, message=FALSE, warning=FALSE}
## data_stage1 wrangling ----
data_stage1$Gender <- as.factor(data_stage1$Gender)
data_stage1$Gender <- plyr::revalue(data_stage1$Gender,
                             c("1" = "Female",
                               "2" = "Male"#,
                               #"3" = "Other"
                               ))
data_stage1$Gender <- factor(data_stage1$Gender,
                         levels = c("Female",
                                    "Male"#,
                                    #"Other"
                                    ))
data_stage1$Gender_SterAct <- as.factor(data_stage1$Gender_SterAct)
data_stage1$Gender_SterAct <- plyr::revalue(data_stage1$Gender_SterAct,
                             c("1" = "Female",
                               "2" = "Male"#,
                               #"3" = "Other"
                               ))
data_stage1$Gender_SterAct <- factor(data_stage1$Gender_SterAct,
                      levels = c("Female",
                                 "Male"#,
                                 #"Other"
                                 ))

data_stage1$gender <- as.factor(ifelse(data_stage1$Gender == "Female" | data_stage1$Gender_SterAct == "Female", 
                      "Female", NA))
#table(data_stage1$gender)
data_stage1$gender <- as.factor(ifelse(is.na(data_stage1$gender) ==TRUE, "Male", "Female"))
#table(data_stage1$gender)
## political interest
# average score
data_stage1$polinterest <- (data_stage1$Pol_Interest_1 + data_stage1$Pol_Interest_2 + data_stage1$Pol_Interest_3 + 
                       data_stage1$Pol_Interest_4 + data_stage1$Pol_Interest_5)/5
## political knowledge
data_stage1$polknowledge <- data_stage1$SC2

# typo in PTK_7 corrected
names(data_stage1)[names(data_stage1) == "PTK_7"] <- "PKT_7"

# recode binary where 1 = correct answer, 0 = incorrect answer or NA
data_stage1$PKT_1_dich <- ifelse(is.na(data_stage1$PKT_1) == TRUE, 0,
         ifelse(data_stage1$PKT_1 == 2, 1, 0))
data_stage1$PKT_2_dich <- ifelse(is.na(data_stage1$PKT_2) == TRUE, 0,
         ifelse(data_stage1$PKT_2 == 1, 1, 0))
data_stage1$PKT_3_dich <- ifelse(is.na(data_stage1$PKT_3) == TRUE, 0,
         ifelse(data_stage1$PKT_3 == 2, 1, 0))
data_stage1$PKT_4_dich <- ifelse(is.na(data_stage1$PKT_4) == TRUE, 0,
         ifelse(data_stage1$PKT_4 == 4, 1, 0))
data_stage1$PKT_5_dich <- ifelse(is.na(data_stage1$PKT_5) == TRUE, 0,
         ifelse(data_stage1$PKT_5 == 1, 1, 0))
data_stage1$PKT_6_dich <- ifelse(is.na(data_stage1$PKT_6) == TRUE, 0,
         ifelse(data_stage1$PKT_6 == 2, 1, 0))
data_stage1$PKT_7_dich <- ifelse(is.na(data_stage1$PKT_7) == TRUE, 0,
         ifelse(data_stage1$PKT_7 == 2, 1, 0))
data_stage1$PKT_8_dich <- ifelse(is.na(data_stage1$PKT_8) == TRUE, 0,
         ifelse(data_stage1$PKT_8 == 1, 1, 0))
data_stage1$PKT_9_dich <- ifelse(is.na(data_stage1$PKT_9) == TRUE, 0,
         ifelse(data_stage1$PKT_9 == 2, 1, 0))
data_stage1$PKT_10_dich <- ifelse(is.na(data_stage1$PKT_10) == TRUE, 0,
         ifelse(data_stage1$PKT_10 == 4, 1, 0))
data_stage1$PKT_11_dich <- ifelse(is.na(data_stage1$PKT_11) == TRUE, 0,
         ifelse(data_stage1$PKT_11 == 1, 1, 0))
data_stage1$PKT_12_dich <- ifelse(is.na(data_stage1$PKT_12) == TRUE, 0,
         ifelse(data_stage1$PKT_12 == 2, 1, 0))
data_stage1$PKT_13_dich <- ifelse(is.na(data_stage1$PKT_13) == TRUE, 0,
         ifelse(data_stage1$PKT_13 == 1, 1, 0)) 
data_stage1$PKT_14_dich <- ifelse(is.na(data_stage1$PKT_14) == TRUE, 0,
         ifelse(data_stage1$PKT_14 == 2, 1, 0))
data_stage1$PKT_15_dich <- ifelse(is.na(data_stage1$PKT_15) == TRUE, 0,
         ifelse(data_stage1$PKT_15 == 3, 1, 0))
data_stage1$PKT_16_dich <- ifelse(is.na(data_stage1$PKT_16) == TRUE, 0,
         ifelse(data_stage1$PKT_16 == 2, 1, 0))
data_stage1$PKT_17_dich <- ifelse(is.na(data_stage1$PKT_17) == TRUE, 0,
         ifelse(data_stage1$PKT_17 == 1, 1, 0))
data_stage1$PKT_18_dich <- ifelse(is.na(data_stage1$PKT_18) == TRUE, 0,
         ifelse(data_stage1$PKT_18 == 1, 1, 0))
data_stage1$PKT_19_dich <- ifelse(is.na(data_stage1$PKT_19) == TRUE, 0,
         ifelse(data_stage1$PKT_19 == 3, 1, 0))
data_stage1$PKT_20_dich <- ifelse(is.na(data_stage1$PKT_20) == TRUE, 0,
         ifelse(data_stage1$PKT_20 == 3, 1, 0))


data_stage1$polknowledge_sum <- 
  data_stage1$PKT_1_dich +
  data_stage1$PKT_2_dich +
  data_stage1$PKT_3_dich +
  data_stage1$PKT_4_dich +
  data_stage1$PKT_5_dich +
  data_stage1$PKT_6_dich +
  data_stage1$PKT_7_dich +
  data_stage1$PKT_8_dich +
  data_stage1$PKT_9_dich +
  data_stage1$PKT_10_dich +
  data_stage1$PKT_11_dich +
  data_stage1$PKT_12_dich +
  data_stage1$PKT_13_dich +
  data_stage1$PKT_14_dich +
  data_stage1$PKT_15_dich +
  data_stage1$PKT_16_dich +
  data_stage1$PKT_17_dich +
  data_stage1$PKT_18_dich +
  data_stage1$PKT_19_dich +
  data_stage1$PKT_20_dich 


data_stage1$diff_PKT <- data_stage1$polknowledge - data_stage1$polknowledge_sum
colsum_diff_PKT <- sum(data_stage1$diff_PKT)
# sum is zero, hence our coding was identical to the automated Qualtrics score
PKT_stage1 <- data_stage1[,307:326]

## condition
data_stage1$condition <- factor(data_stage1$condition,
                         levels = c("control",
                                    "gender_question",
                                    "gender_statement"))
#table(data_stage1$condition)
## politics vs. non-politics


# include in the Politics group individuals who select Politics (=1) or History (=5) or 
# Public Administration (=12) or Government (=27) or Sociology (=11) 
# as topics important for their studies/work

data_stage1$Politics <- ifelse(grepl("1", data_stage1$Field_1) == TRUE  |
                          grepl("1", data_stage1$Field_5) == TRUE  |
                          grepl("1", data_stage1$Field_11) == TRUE  |
                          grepl("1", data_stage1$Field_12) == TRUE  |
                          grepl("1", data_stage1$Field_27) == TRUE, "Yes", "No")

data_stage1$Politics <- as.factor(data_stage1$Politics)
#table(data$Politics)

```

# Results Stage 1

```{r ANCOVA stage 1a}
# Following the analyses reported in the original study and the analysis script made available by the original authors, we tested the replication claim that activation of gender stereotypes influence performance in a political knowledge test with a 2 (gender) × 2 (field of work/study) × 3 (Gender Stereotype Activation) ANCOVA. The dependent variable was participants’ total score on the political knowledge test. As in the original study, a single score of political interest was calculated per participant (i.e., average of responses in the short scale of political interest) and included as a covariate.

ancova <- aov(polknowledge ~ polinterest + gender*Politics*condition, data_stage1)
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

#sjstats::anova_stats(ancovafit)
sjstats::eta_sq(ancovafit, partial = T, ci.lvl = .95)

apa.aov.table(
ancova,
"sample1_h2.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)
```

# Data Stage 2

```{r stage 2 read data, echo=F, message=FALSE, warning=FALSE}
# Data from Stage 2
data_pooled <- sjlabelled::read_spss("analytical_data_stage2.sav")
```

```{r stage 2 exclude cases, echo=F, warning=FALSE}
## Data exclusion ----
data_pooled <- data_pooled[which(data_pooled$Complete == "Yes"),]
```

```{r stage 2 data_pooled wrangling, echo=F, warning=FALSE}
## data_pooled wrangling ----
data_pooled$Gender <- as.factor(data_pooled$Gender)
data_pooled$Gender <- plyr::revalue(data_pooled$Gender,
                             c("1" = "Female",
                               "2" = "Male"#,
                               #"3" = "Other"
                               ))
data_pooled$Gender <- factor(data_pooled$Gender,
                         levels = c("Female",
                                    "Male"#,
                                    #"Other"
                                    ))
data_pooled$Gender_SterAct <- as.factor(data_pooled$Gender_SterAct)
data_pooled$Gender_SterAct <- plyr::revalue(data_pooled$Gender_SterAct,
                             c("1" = "Female",
                               "2" = "Male"#,
                               #"3" = "Other"
                               ))
data_pooled$Gender_SterAct <- factor(data_pooled$Gender_SterAct,
                      levels = c("Female",
                                 "Male"#,
                                 #"Other"
                                 ))

data_pooled$gender <- as.factor(ifelse(data_pooled$Gender == "Female" | data_pooled$Gender_SterAct == "Female", 
                      "Female", NA))
#table(data_pooled$gender)
data_pooled$gender <- as.factor(ifelse(is.na(data_pooled$gender) ==TRUE, "Male", "Female"))
#table(data_pooled$gender)
## political interest
# average score

data_pooled$polinterest <- (data_pooled$Pol_Interest_1 + data_pooled$Pol_Interest_2 + data_pooled$Pol_Interest_3 + 
                       data_pooled$Pol_Interest_4 + data_pooled$Pol_Interest_5)/5
## political knowledge
data_pooled$polknowledge <- data_pooled$SC2

# typo in PTK_7 corrected
names(data_pooled)[names(data_pooled) == "PTK_7"] <- "PKT_7"

# recode binary where 1 = correct answer, 0 = incorrect answer or NA
data_pooled$PKT_1_dich <- ifelse(is.na(data_pooled$PKT_1) == TRUE, 0,
         ifelse(data_pooled$PKT_1 == 2, 1, 0))
data_pooled$PKT_2_dich <- ifelse(is.na(data_pooled$PKT_2) == TRUE, 0,
         ifelse(data_pooled$PKT_2 == 1, 1, 0))
data_pooled$PKT_3_dich <- ifelse(is.na(data_pooled$PKT_3) == TRUE, 0,
         ifelse(data_pooled$PKT_3 == 2, 1, 0))
data_pooled$PKT_4_dich <- ifelse(is.na(data_pooled$PKT_4) == TRUE, 0,
         ifelse(data_pooled$PKT_4 == 4, 1, 0))
data_pooled$PKT_5_dich <- ifelse(is.na(data_pooled$PKT_5) == TRUE, 0,
         ifelse(data_pooled$PKT_5 == 1, 1, 0))
data_pooled$PKT_6_dich <- ifelse(is.na(data_pooled$PKT_6) == TRUE, 0,
         ifelse(data_pooled$PKT_6 == 2, 1, 0))
data_pooled$PKT_7_dich <- ifelse(is.na(data_pooled$PKT_7) == TRUE, 0,
         ifelse(data_pooled$PKT_7 == 2, 1, 0))
data_pooled$PKT_8_dich <- ifelse(is.na(data_pooled$PKT_8) == TRUE, 0,
         ifelse(data_pooled$PKT_8 == 1, 1, 0))
data_pooled$PKT_9_dich <- ifelse(is.na(data_pooled$PKT_9) == TRUE, 0,
         ifelse(data_pooled$PKT_9 == 2, 1, 0))
data_pooled$PKT_10_dich <- ifelse(is.na(data_pooled$PKT_10) == TRUE, 0,
         ifelse(data_pooled$PKT_10 == 4, 1, 0))
data_pooled$PKT_11_dich <- ifelse(is.na(data_pooled$PKT_11) == TRUE, 0,
         ifelse(data_pooled$PKT_11 == 1, 1, 0))
data_pooled$PKT_12_dich <- ifelse(is.na(data_pooled$PKT_12) == TRUE, 0,
         ifelse(data_pooled$PKT_12 == 2, 1, 0))
data_pooled$PKT_13_dich <- ifelse(is.na(data_pooled$PKT_13) == TRUE, 0,
         ifelse(data_pooled$PKT_13 == 1, 1, 0)) 
data_pooled$PKT_14_dich <- ifelse(is.na(data_pooled$PKT_14) == TRUE, 0,
         ifelse(data_pooled$PKT_14 == 2, 1, 0))
data_pooled$PKT_15_dich <- ifelse(is.na(data_pooled$PKT_15) == TRUE, 0,
         ifelse(data_pooled$PKT_15 == 3, 1, 0))
data_pooled$PKT_16_dich <- ifelse(is.na(data_pooled$PKT_16) == TRUE, 0,
         ifelse(data_pooled$PKT_16 == 2, 1, 0))
data_pooled$PKT_17_dich <- ifelse(is.na(data_pooled$PKT_17) == TRUE, 0,
         ifelse(data_pooled$PKT_17 == 1, 1, 0))
data_pooled$PKT_18_dich <- ifelse(is.na(data_pooled$PKT_18) == TRUE, 0,
         ifelse(data_pooled$PKT_18 == 1, 1, 0))
data_pooled$PKT_19_dich <- ifelse(is.na(data_pooled$PKT_19) == TRUE, 0,
         ifelse(data_pooled$PKT_19 == 3, 1, 0))
data_pooled$PKT_20_dich <- ifelse(is.na(data_pooled$PKT_20) == TRUE, 0,
         ifelse(data_pooled$PKT_20 == 3, 1, 0))


data_pooled$polknowledge_sum <- 
  data_pooled$PKT_1_dich +
  data_pooled$PKT_2_dich +
  data_pooled$PKT_3_dich +
  data_pooled$PKT_4_dich +
  data_pooled$PKT_5_dich +
  data_pooled$PKT_6_dich +
  data_pooled$PKT_7_dich +
  data_pooled$PKT_8_dich +
  data_pooled$PKT_9_dich +
  data_pooled$PKT_10_dich +
  data_pooled$PKT_11_dich +
  data_pooled$PKT_12_dich +
  data_pooled$PKT_13_dich +
  data_pooled$PKT_14_dich +
  data_pooled$PKT_15_dich +
  data_pooled$PKT_16_dich +
  data_pooled$PKT_17_dich +
  data_pooled$PKT_18_dich +
  data_pooled$PKT_19_dich +
  data_pooled$PKT_20_dich 


data_pooled$diff_PKT <- data_pooled$polknowledge - data_pooled$polknowledge_sum
sum(data_pooled$diff_PKT)
# sum is zero, hence our coding was identical to the automated Qualtrics score
PKT_pooled <- data_pooled[,307:326]
## condition
data_pooled$condition <- factor(data_pooled$condition,
                         levels = c("control",
                                    "gender_question",
                                    "gender_statement"))
#table(data_pooled$condition)
## politics vs. non-politics
data_pooled$Politics <- factor(data_pooled$Politics,
                         levels = c("Yes",
                                    "No"))

# include in the Politics group individuals who select Politics (=1) or History (=5) or 
# Public Administration (=12) or Government (=27) or Sociology (=11) 
# as topics important for their studies/work

data_pooled$Politics <- ifelse(grepl("1", data_pooled$Field_1) == TRUE  |
                          grepl("1", data_pooled$Field_5) == TRUE  |
                          grepl("1", data_pooled$Field_11) == TRUE  |
                          grepl("1", data_pooled$Field_12) == TRUE  |
                          grepl("1", data_pooled$Field_27) == TRUE, "Yes", "No")
#table(data_pooled$Politics)

```

# Results Stage 2

The pooled analytical sample (first and second stages of data collection together) was composed by *N* = `r nrow(data_pooled)` participants (*M age* = `r round(mean(data_pooled$Age),2)` years, *SD* = `r round(sd(data_pooled$Age),2)`, `r round((nrow(data_pooled[which(data_pooled$gender == "Female"),])/nrow(data_pooled))*100,2)`% female). 

```{r stage 2 sample composition A, echo=F, message=FALSE, warning=FALSE}
psych::describe(data_pooled$Age)
round((nrow(data_pooled[which(data_pooled$gender == "Female"),])/nrow(data_pooled))*100,2)
```

```{r ANCOVA stage 2a}
# 2 (gender) × 2 (field of work/study) × 3 (Gender Stereotype Activation) ANCOVA. The dependent variable was participants’ total score on the political knowledge test. As in the original study, a single score of political interest was calculated per participant (i.e., average of responses in the short scale of political interest) and included as a covariate.

ancova <- aov(polknowledge ~ polinterest + gender*Politics*condition, data_pooled)
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

psych::describeBy(data_pooled$polknowledge, group = data_pooled$gender)

#sjstats::anova_stats(ancovafit)
sjstats::eta_sq(ancovafit, partial = T, ci.lvl = .95)

apa.aov.table(
ancova,
"sample2_h2.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)
```

```{r emmeans stage 2a}
# Bonferroni-corrected pairwise comparisons of gender with the emmeans function in R (Lenth et al. 2018).
data_pooled$polinterest_c <- center_scale(data_pooled$polinterest)

ancova <- aov(polknowledge ~ polinterest_c + gender*Politics*condition, data_pooled)
emmeans.ancova <- emmeans(ancova, specs = c("gender", "condition"))
emmeans.ancova 

contrast <- contrast(emmeans.ancova, "pairwise", 
                        simple = "gender", # "each" gives all contrasts
                        combine = F, 
                        adjust = "bonferroni")

contrast
```

```{r plot gender differences stage 2}

symnum.args <- list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 1), symbols = c("****", "***", "**", "*", "ns"))

# Visualization of gender differences in political knowledge test per stereotype activation condition
p <- ggplot(data_pooled, aes(x = gender, y = polknowledge, 
                             color = gender))  +
  scale_color_manual(values = c("#66a182", "#d1495b")) +
  #scale_color_manual(values = c("black", "grey40")) + # b/w version
  geom_boxplot() +
  facet_wrap(~condition) +
  labs(
    y = "Political Knowledge\n",
    x = ""
    #title = "Mean Comparison of Political Knowledge Test Between Males and Females per \nExperimental Condition.\n"
  ) +
  geom_jitter(shape=16, 
              position=position_jitter(0.2),
              alpha = .4) +
  theme_bw() +
  stat_compare_means(method = "t.test", 
                       label = "p.signif",
                       label.x = 1.5,
                       symnum.args = symnum.args,
                       family = "serif",
                     size = 5) + 
     stat_summary(fun.data = function(x) data.frame(y=-1, label = paste("M =",round(mean(x),2))),
                  geom="text",
                  inherit.aes = T,
                  color = "black",
                  family = "serif",
                  size = 5) +
     theme(legend.position="none", 
           strip.text = element_text(size = 15),
           text = element_text(size=15, family="serif"),
           axis.text.x = element_text(color = "black", size = 15))

p
```

# Exploratory Analyses

## Evidence-updated replication Bayes factors
Evidence-updated replication Bayes factors were calculated in JASP.

## Scale difficulty z-test

```{r z-test}
# z-test Stage 1
x1 <- mean(data_stage1$polknowledge)
x2 <- mean(originaldata$PolWis_score)

sigma1 <- sd(data_stage1$polknowledge)/sqrt(nrow(data_stage1))
sigma2 <- sd(originaldata$PolWis_score)/sqrt(nrow(originaldata))

z.test(x1, x2, sigma1, sigma2)

# z-Test Stage 2
x1 <- mean(data_pooled$polknowledge)
x2 <- mean(originaldata$PolWis_score)

sigma1 <- sd(data_pooled$polknowledge)/sqrt(nrow(data_pooled))
sigma2 <- sd(originaldata$PolWis_score)/sqrt(nrow(originaldata))

z.test(x1, x2, sigma1, sigma2)

# plot political knowledge score distributions
p <- ggplot(data_stage1, aes(x = polknowledge)) + 
    geom_density(alpha=.5, fill="#00798c") +
  theme_bw() +
  labs(
    y = "Density",
    x = "Political Knowledge Score",
    title = "Comparison of Political Knowledge Test Density Plots of the Original and the Replication Data.",
    subtitle = "Stage 1 (Upper Panel) and Stage 2 (Lower Panel)\n"
  ) +
  theme(text = element_text(family = "serif", color = "black"))

p <- p + geom_density(mapping = aes(x = PolWis_score),
                 data = originaldata,
                 alpha = .5,
                 fill = "#d1495b") +
  annotate("text", x = 13, y = 0.09, label = "Replication Data Stage 1", family = "serif") +
  annotate("text", x = 17, y = 0.05, label = "Original Data", family = "serif")

#p

p2 <- ggplot(data_pooled, aes(x = polknowledge)) + 
    geom_density(alpha=.5, fill="#00798c") +
  theme_bw() +
  labs(
    y = "Density",
    x = "Political Knowledge Score",
    caption = "Original Data in Red, Replication Data in Blue"
  ) +
  theme(text = element_text(family = "serif", color = "black"))

p2 <- p2 + geom_density(mapping = aes(x = PolWis_score),
                 data = originaldata,
                 alpha = .5,
                 fill = "#d1495b")+
  annotate("text", x = 13, y = 0.09, label = "Replication Data Stage 2", family = "serif") +
  annotate("text", x = 17, y = 0.05, label = "Original Data", family = "serif")

#p2

grid.arrange(p, p2)
```

## Scale difficulty IRT

```{r IRT 2PL}
# Original data

# recode binary where 1 = correct answer, 0 = incorrect answer or NA
originaldata$PKT_1_dich <- ifelse(is.na(originaldata$polwis01) == TRUE, 0,
                                 ifelse(originaldata$polwis01 == 4, 1, 0))
originaldata$PKT_2_dich <- ifelse(is.na(originaldata$polwis02) == TRUE, 0,
                                 ifelse(originaldata$polwis02 == 4, 1, 0))
originaldata$PKT_3_dich <- ifelse(is.na(originaldata$polwis03) == TRUE, 0,
                                 ifelse(originaldata$polwis03 == 1, 1, 0))
originaldata$PKT_4_dich <- ifelse(is.na(originaldata$polwis04) == TRUE, 0,
                                 ifelse(originaldata$polwis04 == 2, 1, 0))
originaldata$PKT_5_dich <- ifelse(is.na(originaldata$polwis05) == TRUE, 0,
                                 ifelse(originaldata$polwis05 == 2, 1, 0))
originaldata$PKT_6_dich <- ifelse(is.na(originaldata$polwis06) == TRUE, 0,
                                 ifelse(originaldata$polwis06 == 3, 1, 0))
originaldata$PKT_7_dich <- ifelse(is.na(originaldata$polwis07) == TRUE, 0,
                                 ifelse(originaldata$polwis07 == 2, 1, 0))
originaldata$PKT_8_dich <- ifelse(is.na(originaldata$polwis08) == TRUE, 0,
                                 ifelse(originaldata$polwis08 == 2, 1, 0))
originaldata$PKT_9_dich <- ifelse(is.na(originaldata$polwis09) == TRUE, 0,
                                 ifelse(originaldata$polwis09 == 4, 1, 0))
originaldata$PKT_10_dich <- ifelse(is.na(originaldata$polwis10) == TRUE, 0,
                                  ifelse(originaldata$polwis10 == 1, 1, 0))
originaldata$PKT_11_dich <- ifelse(is.na(originaldata$polwis11) == TRUE, 0,
                                  ifelse(originaldata$polwis11 == 1, 1, 0))
originaldata$PKT_12_dich <- ifelse(is.na(originaldata$polwis12) == TRUE, 0,
                                  ifelse(originaldata$polwis12 == 3, 1, 0))
originaldata$PKT_13_dich <- ifelse(is.na(originaldata$polwis13) == TRUE, 0,
                                  ifelse(originaldata$polwis13 == 3, 1, 0)) 
originaldata$PKT_14_dich <- ifelse(is.na(originaldata$polwis14) == TRUE, 0,
                                  ifelse(originaldata$polwis14 == 2, 1, 0))
originaldata$PKT_15_dich <- ifelse(is.na(originaldata$polwis15) == TRUE, 0,
                                  ifelse(originaldata$polwis15 == 1, 1, 0))
originaldata$PKT_16_dich <- ifelse(is.na(originaldata$polwis16) == TRUE, 0,
                                  ifelse(originaldata$polwis16 == 1, 1, 0))
originaldata$PKT_17_dich <- ifelse(is.na(originaldata$polwis17) == TRUE, 0,
                                  ifelse(originaldata$polwis17 == 4, 1, 0))
originaldata$PKT_18_dich <- ifelse(is.na(originaldata$polwis18) == TRUE, 0,
                                  ifelse(originaldata$polwis18 == 2, 1, 0))
originaldata$PKT_19_dich <- ifelse(is.na(originaldata$polwis19) == TRUE, 0,
                                  ifelse(originaldata$polwis19 == 3, 1, 0))
originaldata$PKT_20_dich <- ifelse(is.na(originaldata$polwis20) == TRUE, 0,
                                  ifelse(originaldata$polwis20 == 2, 1, 0))

#
PolKnow_items_Orig <- originaldata[, c("PKT_1_dich", "PKT_2_dich", 
                                       "PKT_3_dich", "PKT_4_dich", "PKT_5_dich", "PKT_6_dich", "PKT_7_dich", 
                                       "PKT_8_dich", "PKT_9_dich", "PKT_10_dich", "PKT_11_dich", "PKT_12_dich", 
                                       "PKT_13_dich", "PKT_14_dich", "PKT_15_dich", "PKT_16_dich", "PKT_17_dich", 
                                       "PKT_18_dich", "PKT_19_dich", "PKT_20_dich")]
#
names(PolKnow_items_Orig) <- c("item 1", "item 2", "item 3", "item 4", "item 5", "item 6", 
                              "item 7", "item 8", "item 9", "item 10", "item 11", "item 12", 
                              "item 13", "item 14", "item 15", "item 16", "item 17", "item 18", 
                              "item 19", "item 20")

######################################
# Pooled data (Stage 1 and 2 together)
######################################

PolKnow_items_all <- data_pooled[, c("PKT_1_dich", "PKT_2_dich", "PKT_3_dich", 
                                    "PKT_4_dich", "PKT_5_dich", "PKT_6_dich",
                                    "PKT_7_dich", "PKT_8_dich", "PKT_9_dich", 
                                    "PKT_10_dich", "PKT_11_dich", "PKT_12_dich", 
                                    "PKT_13_dich", "PKT_14_dich",
                                    "PKT_15_dich", "PKT_16_dich", "PKT_17_dich", 
                                    "PKT_18_dich", "PKT_19_dich", "PKT_20_dich")]
#
names(PolKnow_items_all) <- c("item 1", "item 2", "item 3", "item 4", "item 5", "item 6", 
                              "item 7", "item 8", "item 9", "item 10", "item 11", "item 12", 
                              "item 13", "item 14", "item 15", "item 16", "item 17", "item 18", 
                              "item 19", "item 20")
#
#
#
# ggcorrplot::ggcorrplot(cor(PolKnow_items_all,  use = "pairwise.complete.obs"), 
#                        hc.order = TRUE, type = "lower", lab = TRUE, insig = "blank")
# #
# ggcorrplot::ggcorrplot(cor(PolKnow_items_Orig,  use = "pairwise.complete.obs"), 
#                        hc.order = TRUE, type = "lower", lab = TRUE, insig = "blank")
#
#
S1.Rasch.OD <- mirt::mirt(PolKnow_items_Orig, 1, itemtype = 'Rasch', SE = TRUE, technical = list(removeEmptyRows=TRUE))
S1.2PL.OD   <- mirt::mirt(PolKnow_items_Orig, 1, itemtype = '2PL', SE = TRUE, technical = list(removeEmptyRows=TRUE))
S1.3PL.OD   <- mirt::mirt(PolKnow_items_Orig, 1, itemtype = '3PL', SE = TRUE, technical = list(removeEmptyRows=TRUE))
#
anova(S1.Rasch.OD, S1.2PL.OD)# 2PL 
anova(S1.2PL.OD, S1.3PL.OD)  # 2PL
#
#
S1.Rasch.pooled <- mirt::mirt(PolKnow_items_all, 1, itemtype = 'Rasch', SE = TRUE, technical = list(removeEmptyRows=TRUE))
S1.2PL.pooled <- mirt::mirt(PolKnow_items_all, 1, itemtype = '2PL', SE = TRUE, technical = list(removeEmptyRows=TRUE))
S1.3PL.pooled <- mirt::mirt(PolKnow_items_all, 1, itemtype = '3PL', SE = TRUE, technical = list(removeEmptyRows=TRUE))
#
anova(S1.Rasch.pooled,S1.2PL.pooled)# 2PL 
anova(S1.2PL.pooled,S1.3PL.pooled)  # 3PL 
#
#
library(gridExtra)
library(latticeExtra)
#
#
#mirt::?plot
#
par(family="Times")
windowsFonts(A = windowsFont("Times New Roman"))
key=list(columns=1, text=list(lab=c("Ihme and Tausendpfund (2018)","Replication")), 
         lines=list(lwd=1, col=c("blue","orange")))
#
#
# Reliability rxx
# -------------------
plot7b <- plot(S1.2PL.OD, type = "rxx", facet_items=F, theta_lim = c(-6, 6), key=key, main="Reliability", family="A")
plot8b <- update(plot(S1.2PL.pooled, type = "rxx", facet_items=F, theta_lim = c(-6, 6)), col="orange", family="A")
plot7b + plot8b
#
#
#
# Expected Total Score  (trace)
# ------------------------------
plot9b <- plot(S1.2PL.OD, type = "score", facet_items=F, theta_lim = c(-6, 6), key=key, main="Expected Total Score")
plot10b <- update(plot(S1.2PL.pooled, type = "score", facet_items=F, theta_lim = c(-6, 6)), col="orange")
plot9b + plot10b
#
#
#
# Item Trace Lines  (faceted)
# ------------------------------
plot11b <- plot(S1.2PL.OD, type = "trace", facet_items=T, theta_lim = c(-6, 6), key=key, main="Item Trace Lines")
plot12b <- update(plot(S1.2PL.pooled, type = "trace", facet_items=T, theta_lim = c(-6, 6)), col="orange")
plot11b + plot12b
#
#
#
library(curl)
library(ggplot2)
library(grid)
library(gridExtra)
#
lay = rbind(c(1,3,3),
            c(2,3,3))

print(grid.arrange(arrangeGrob(plot7b + plot8b, left = textGrob("a)", x = unit(1, "npc"), 
                                                       y = unit(.95, "npc"))), 
                   arrangeGrob(plot9b + plot10b, left =textGrob("b)", x = unit(1, "npc"), 
                                                      y = unit(1, "npc"))),
                   arrangeGrob(plot11b + plot12b, left=textGrob("c)", x = unit(1, "npc"), 
                                                     y = unit(.95, "npc"))),
                   layout_matrix = lay))
```

# Supplementary Materials

## Reliability political knowledge scale

```{r}
cronbach.alpha(PKT_stage1, CI = T)
cronbach.alpha(PKT_pooled, CI = T)
```

## First stage of data collection - detailed results
### Sample Stage 1

```{r stage 1 sample composition, echo=F, message=FALSE, warning=FALSE}
# read in raw data
data_stage1_raw <- sjlabelled::read_spss("analytical_data_stage1.sav")

# participants who did not consent (Yes = 1, No = 2)
table(data_stage1_raw$Consent, useNA = "ifany")
data_stage1_raw <- data_stage1_raw[which(data_stage1_raw$Consent == "1"),]

# participants who did not commit to answering truthfully
table(data_stage1_raw$Data_Quality, useNA = "ifany")
data_stage1_raw <- data_stage1_raw[which(data_stage1_raw$Data_Quality == "1"),]

# participants who failed the quota
table(data_stage1_raw$Failed_Quota, useNA = "ifany")
data_stage1_raw <- data_stage1_raw[which(data_stage1_raw$Failed_Quota == "No"),]

# number of participants allocated to conditions
table(data_stage1_raw$condition, useNA = "ifany")

# participants who indicated that they had cheated or used help on the political knowledge test
table(data_stage1_raw$condition, data_stage1_raw$Help_Ihme, useNA = "ifany")
data_stage1_raw <- data_stage1_raw[which(data_stage1_raw$Help_Ihme == "1"),]

# participants who failed attention check questions
table(data_stage1_raw$condition, data_stage1_raw$Attention_Score, useNA = "ifany")
data_stage1_raw <- data_stage1_raw[which(data_stage1_raw$Attention_Score == "0"),]

# participants who identified as neither male nor female
table(data_stage1_raw$condition,data_stage1_raw$Gender_SterAct, useNA = "ifany")
table(data_stage1_raw$condition,data_stage1_raw$Gender, useNA = "ifany")

# participants that are younger than 18
table(data_stage1_raw$condition,data_stage1_raw$Age, useNA = "ifany")
data_stage1_raw <- data_stage1_raw[-which(data_stage1_raw$Age < 18),]

table(data_stage1_raw$condition, useNA = "ifany")
```

### Frequency of participants per conditions

Table S2
```{r stage 1 crosstables per condition, echo=F, message=FALSE, warning=FALSE}
prop.table(table(data_stage1$Politics, useNA = "ifany"))
prop.table(table(data_stage1$gender, useNA = "ifany"))
prop.table(table(data_stage1$gender, data_stage1$condition, data_stage1$Politics, useNA = "ifany"))
```

### Statistical analyses

Table S3
```{r ANOVA stage 1b}
# 2 (gender: male vs female) × 2 (field of work/study: Political vs Non-Political) × 3 (Gender Stereotype Activation: stereotype not activated vs. stereotype activated by gender-question vs. stereotype activated by gender difference statement) ANOVA with participants’ score in the political knowledge test as the dependent variable. 

aov <- aov(polknowledge ~ gender*Politics*condition, data_stage1)
aovfit <- car::Anova(aov, digits = 3, type = 3)
aovfit

#sjstats::anova_stats(aovfit, digits = 3, partial = T)
sjstats::eta_sq(aovfit, partial = T, ci.lvl = .95)

apa.aov.table(
aov,
"sample1_h1.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)
```

Table S4
```{r ANCOVA stage 1}
# 2 (gender) × 2 (field of work/study) × 3 (Gender Stereotype Activation) ANCOVA. The dependent variable was participants’ total score on the political knowledge test. As in the original study, a single score of political interest was calculated per participant (i.e., average of responses in the short scale of political interest) and included as a covariate.

ancova <- aov(polknowledge ~ polinterest + gender*Politics*condition, data_stage1)
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

#sjstats::anova_stats(ancovafit)
sjstats::eta_sq(ancovafit, partial = T, ci.lvl = .95)

apa.aov.table(
ancova,
"sample1_h2.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)
```

```{r emmeans stage 1}
# Bonferroni-corrected pairwise comparisons of gender X Politics interaction with the emmeans function in R (Lenth et al. 2018).
ancova <- aov(polknowledge ~ polinterest + gender*Politics*condition, data_stage1)
emmeans.ancova <- emmeans(ancova, specs = c("gender", "Politics"))
emmeans.ancova 

contrast <- contrast(emmeans.ancova, "pairwise", 
                        simple = "Politics", # "each" gives all contrasts
                        combine = F, 
                        adjust = "bonferroni")

contrast

# Bonferroni-corrected pairwise comparisons of Politics X condition interaction with the emmeans function in R (Lenth et al. 2018).
ancova <- aov(polknowledge ~ polinterest + gender*Politics*condition, data_stage1)
emmeans.ancova <- emmeans(ancova, specs = c("condition", "Politics"))
emmeans.ancova 

contrast <- contrast(emmeans.ancova, "pairwise", 
                        simple = "Politics", # "each" gives all contrasts
                        combine = F, 
                        adjust = "bonferroni")

contrast

# Bonferroni-corrected pairwise comparisons of gender X Politics x condition interaction with the emmeans function in R (Lenth et al. 2018).
ancova <- aov(polknowledge ~ polinterest + gender*Politics*condition, data_stage1)
emmeans.ancova <- emmeans(ancova, specs = c("gender", "Politics", "condition"))
emmeans.ancova 

contrast <- contrast(emmeans.ancova, "pairwise", 
                        simple = "gender", # "each" gives all contrasts
                        combine = F, 
                        adjust = "bonferroni")

contrast
```

## Detailed results for the pooled data

### Sample Stage 2

```{r stage 2 sample composition, echo=F, message=FALSE, warning=FALSE}
# read in raw data
data_stage2_raw <- read.csv("data_stage2.csv")

# participants who did not consent (Yes = 1, No = 2)
table(data_stage2_raw$Consent, useNA = "ifany")
data_stage2_raw <- data_stage2_raw[which(data_stage2_raw$Consent == "1"),]

# participants who did not commit to answering truthfully
table(data_stage2_raw$Data_Quality, useNA = "ifany")
data_stage2_raw <- data_stage2_raw[which(data_stage2_raw$Data_Quality == "1"),]

# participants who failed the quota
table(data_stage2_raw$Failed_Quota, useNA = "ifany")
data_stage2_raw <- data_stage2_raw[which(data_stage2_raw$Failed_Quota == "No"),]

# number of participants allocated to conditions
table(data_stage2_raw$condition, useNA = "ifany")

# participants who didn't finish their questionnaire
table(data_stage2_raw$Finished,data_stage2_raw$condition, useNA = "ifany")
data_stage2_raw <- data_stage2_raw[which(data_stage2_raw$Finished == "1"),]

# participants who indicated that they had cheated or used help on the political knowledge test
table(data_stage2_raw$condition, data_stage2_raw$Help_Ihme, useNA = "ifany")
data_stage2_raw <- data_stage2_raw[which(data_stage2_raw$Help_Ihme == "1"),]

# participants who failed attention check questions
table(data_stage2_raw$condition, data_stage2_raw$Attention_Score, useNA = "ifany")
data_stage2_raw <- data_stage2_raw[which(data_stage2_raw$Attention_Score == "0"),]

# participants who identified as neither male nor female
table(data_stage2_raw$condition,data_stage2_raw$Gender_SterAct, useNA = "ifany")
table(data_stage2_raw$condition,data_stage2_raw$Gender, useNA = "ifany")
data_stage2_raw <- data_stage2_raw[-which(data_stage2_raw$Gender_SterAct == "3"),]
data_stage2_raw <- data_stage2_raw[-which(data_stage2_raw$Gender == "3"),]

# participants that are younger than 18
table(data_stage2_raw$condition,data_stage2_raw$Age, useNA = "ifany")
data_stage2_raw <- data_stage2_raw[-which(data_stage2_raw$Age == 6),]

table(data_stage2_raw$condition, useNA = "ifany")
```

### Frequency of participants per conditions

Table S5
```{r stage 2 crosstables per condition, echo=F, message=FALSE, warning=FALSE}
prop.table(table(data_pooled$Politics, useNA = "ifany"))
prop.table(table(data_pooled$gender, useNA = "ifany"))
prop.table(table(data_pooled$gender, data_pooled$condition, data_pooled$Politics, useNA = "ifany"))
```

### Statistical analyses

Table S6
```{r ANOVA stage 2}
# 2 (gender: male vs female) × 2 (field of work/study: Political vs Non-Political) × 3 (Gender Stereotype Activation: stereotype not activated vs. stereotype activated by gender-question vs. stereotype activated by gender difference statement) ANOVA with participants’ score in the political knowledge test as the dependent variable. 

aov <- aov(polknowledge ~ gender*Politics*condition, data_pooled)
aovfit <- car::Anova(aov, digits = 3, type = 3)
aovfit

#sjstats::anova_stats(aovfit, digits = 3, partial = T)
sjstats::eta_sq(aovfit, partial = T, ci.lvl = .95)

apa.aov.table(
aov,
"sample2_h1.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)
```

Table S7
```{r ANCOVA stage 2b}
# 2 (gender) × 2 (field of work/study) × 3 (Gender Stereotype Activation) ANCOVA. The dependent variable was participants’ total score on the political knowledge test. As in the original study, a single score of political interest was calculated per participant (i.e., average of responses in the short scale of political interest) and included as a covariate.

ancova <- aov(polknowledge ~ polinterest + gender*Politics*condition, data_pooled)
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

#sjstats::anova_stats(ancovafit)
sjstats::eta_sq(ancovafit, partial = T, ci.lvl = .95)

apa.aov.table(
ancova,
"sample2_h2.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)

# 2 x 2 x 3 ANCOVA with original data

ancova <- aov(PolWis_score ~ PolInt_score + gender_alle*field_of_study*Stereotypaktivierung, originaldata)
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

#sjstats::anova_stats(ancovafit)
sjstats::eta_sq(ancovafit, partial = T, ci.lvl = .95)

apa.aov.table(
ancova,
"original_h2.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)


```

```{r emmeans stage 2}
# Bonferroni-corrected pairwise comparisons of gender X Politics interaction with the emmeans function in R (Lenth et al. 2018).
ancova <- aov(polknowledge ~ polinterest + gender*Politics*condition, data_pooled)
emmeans.ancova <- emmeans(ancova, specs = c("gender", "condition"))
emmeans.ancova 

contrast <- contrast(emmeans.ancova, "pairwise", 
                        simple = "condition", # "each" gives all contrasts
                        combine = F, 
                        adjust = "bonferroni")

contrast
```

```{r one-tailed t-tests stage 2}
## One-tailed t-tests like in the original paper
# WOMEN
# control - question
t.test(data_pooled[which(data_pooled$condition != "gender_question" & data_pooled$gender == "Female"),]$polknowledge ~ data_pooled[which(data_pooled$condition != "gender_question" & data_pooled$gender == "Female"),]$condition,
       alternative = "less")
#
effectsize::cohens_d(data_pooled[which(data_pooled$condition != "gender_question" & data_pooled$gender == "Female"),]$polknowledge ~ data_pooled[which(data_pooled$condition != "gender_question" & data_pooled$gender == "Female"),]$condition,
       alternative = "less")
# control - statement
t.test(data_pooled[which(data_pooled$condition != "gender_statement" & data_pooled$gender == "Female"),]$polknowledge ~ data_pooled[which(data_pooled$condition != "gender_statement" & data_pooled$gender == "Female"),]$condition,
       alternative = "less")
#
effectsize::cohens_d(data_pooled[which(data_pooled$condition != "gender_statement" & data_pooled$gender == "Female"),]$polknowledge ~ data_pooled[which(data_pooled$condition != "gender_statement" & data_pooled$gender == "Female"),]$condition,
       alternative = "less")
# MEN
# control - question
t.test(data_pooled[which(data_pooled$condition != "gender_question" & data_pooled$gender == "Male"),]$polknowledge ~ data_pooled[which(data_pooled$condition != "gender_question" & data_pooled$gender == "Male"),]$condition,
       alternative = "greater")
#
effectsize::cohens_d(data_pooled[which(data_pooled$condition != "gender_question" & data_pooled$gender == "Male"),]$polknowledge ~ data_pooled[which(data_pooled$condition != "gender_question" & data_pooled$gender == "Male"),]$condition,
       alternative = "greater")
# control - statement
t.test(data_pooled[which(data_pooled$condition != "gender_statement" & data_pooled$gender == "Male"),]$polknowledge ~ data_pooled[which(data_pooled$condition != "gender_statement" & data_pooled$gender == "Male"),]$condition,
       alternative = "greater")
#
effectsize::cohens_d(data_pooled[which(data_pooled$condition != "gender_statement" & data_pooled$gender == "Male"),]$polknowledge ~ data_pooled[which(data_pooled$condition != "gender_statement" & data_pooled$gender == "Male"),]$condition,
       alternative = "greater")
```

## Exploratory Analyses
### Analyses of “don’t know” and incorrect answers in the political knowledge test

Table S8 and S9
```{r Dont know answers in the political knowledge test ANCOVA stage 2}
#Frequency of don't know answers per participant
# In qualtrics, "don't know responses are coded as 0"

dput(names(data_pooled))
data_pooled$PKT_DK <- rowSums(data_pooled[, c("PKT_1", "PKT_2", "PKT_3", "PKT_4", "PKT_5", "PKT_6", "PKT_7", "PKT_8", "PKT_9", "PKT_10", "PKT_11", "PKT_12", "PKT_13", "PKT_14", "PKT_15", "PKT_16", "PKT_17", "PKT_18", "PKT_19", "PKT_20" )] == 0, na.rm = T)

ancova <- aov(PKT_DK ~ polinterest + gender*condition, data_pooled)
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

sjstats::anova_stats(ancovafit)
sjstats::eta_sq(ancovafit, partial = T, ci.lvl = .95)

#Table of means and sd of don't know responses per gender, per condition

group_by(data_pooled, gender, condition) %>%
  summarise(
    count = n(),
    mean = mean(PKT_DK, na.rm = TRUE),
    sd = sd(PKT_DK, na.rm = TRUE)
  )

#Table of means and sd of don't know responses per gender

group_by(data_pooled, gender)%>%
  summarise(
    count = n(),
    mean = mean(PKT_DK, na.rm = TRUE),
    sd = sd(PKT_DK, na.rm = TRUE)
  )
#
#
apa.aov.table(
ancova,
"sample2_ex_DK.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)

```

Table S10 and S11
```{r Incorrect answers in the political knowledge test stage 2}
data_pooled$PKT_incorrect <- rowSums(data_pooled[, c("PKT_1", "PKT_2", "PKT_3", "PKT_4", "PKT_5", "PKT_6", "PKT_7", "PKT_8", "PKT_9", 
                            "PKT_10", "PKT_11", "PKT_12", "PKT_13", "PKT_14", "PKT_15", "PKT_16", "PKT_17", 
                            "PKT_18", "PKT_19", "PKT_20" )] == 2, na.rm = T)
#Frequency of don't know answers per participant
# In qualtrics, "don't know responses are coded as 0"
ancova <- aov(PKT_incorrect ~ polinterest + gender*condition, data_pooled)
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

sjstats::anova_stats(ancovafit)
sjstats::eta_sq(ancovafit, partial = T, ci.lvl = .95)

#Table of means and sd of incorrect responses per gender, per condition

group_by(data_pooled, condition, gender) %>%
  summarise(
    count = n(),
    mean = mean(PKT_incorrect, na.rm = TRUE),
    sd = sd(PKT_incorrect, na.rm = TRUE)
  )

#Table of means and sd of incorrect responses per gender
group_by(data_pooled, gender)%>%
  summarise(
    count = n(),
    mean = mean(PKT_incorrect, na.rm = TRUE),
    sd = sd(PKT_incorrect, na.rm = TRUE)
  )
#

emmeans.ancova <- emmeans(ancova, specs = c("gender", "condition"))
emmeans.ancova 

contrast <- contrast(emmeans.ancova, "pairwise", 
                        simple = "gender", # "each" gives all contrasts
                        combine = F, 
                        adjust = "bonferroni")

contrast

#
apa.aov.table(
ancova,
"sample2_ex_inc.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)
```

Table S12
```{r DK/incorrect proportion as covariate in the political knowledge test ANCOVA stage 2}
data_pooled$prop.dk <- data_pooled$PKT_DK/(data_pooled$PKT_DK + data_pooled$PKT_incorrect)

ancova <- aov(polknowledge ~ polinterest + prop.dk + gender*Politics*condition, data_pooled)
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

sjstats::anova_stats(ancovafit)
sjstats::eta_sq(ancovafit, partial = T, ci.lvl = .95)

#Table of means and sd of incorrect responses per gender, per condition

group_by(data_pooled, condition, gender) %>%
  summarise(
    count = n(),
    mean = mean(PKT_DK, na.rm = TRUE),
    sd = sd(PKT_DK, na.rm = TRUE)
  )

#Table of means and sd of incorrect responses per gender
group_by(data_pooled, gender)%>%
  summarise(
    count = n(),
    mean = mean(PKT_DK, na.rm = TRUE),
    sd = sd(PKT_DK, na.rm = TRUE)
  )
#
emmeans.ancova <- emmeans(ancova, specs = c("gender", "condition"))
emmeans.ancova 

contrast <- contrast(emmeans.ancova, "pairwise", 
                        simple = "gender", # "each" gives all contrasts
                        combine = F, 
                        adjust = "bonferroni")

contrast

#
apa.aov.table(
ancova,
"sample2_ex_prop.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)
```

### Analyses of attempted answers in the political knowledge test

Table S13 and S14
```{r Attempted items in the political knowledge test stage 2}

#Frequency of attempted items per participant (items that were not answered with "don't know" and were not skipped)

data_pooled$attempt <- rowSums(data_pooled[, c("PKT_1", "PKT_2", "PKT_3", "PKT_4", "PKT_5", "PKT_6", "PKT_7", "PKT_8", "PKT_9", 
                            "PKT_10", "PKT_11", "PKT_12", "PKT_13", "PKT_14", "PKT_15", "PKT_16", "PKT_17", 
                            "PKT_18", "PKT_19", "PKT_20" )] != 0, na.rm = T)

ancova <- aov(attempt ~ polinterest + gender*condition, data_pooled)
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

sjstats::anova_stats(ancovafit)
sjstats::eta_sq(ancovafit, partial = T, ci.lvl = .95)

#Table of means and sd of attempted responses per gender, per condition

group_by(data_pooled, gender, condition) %>%
  summarise(
    count = n(),
    mean = round(mean(attempt, na.rm = TRUE), digits = 3),
    sd = sd(attempt, na.rm = TRUE)
  )

#Table of means and sd of attempted responses per gender

group_by(data_pooled, gender)%>%
  summarise(
    count = n(),
    mean = mean(attempt, na.rm = TRUE),
    sd = sd(attempt, na.rm = TRUE))
#
apa.aov.table(
ancova,
"sample2_ex_attempt.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)
```

### Political interest

Table S15
```{r Political interest stage 2}
aov <- aov(polinterest ~ gender*Politics*condition, data_pooled)
aovfit <- car::Anova(aov, type = 3)
aovfit

sjstats::anova_stats(aovfit)
sjstats::eta_sq(aovfit, partial = T, ci.lvl = .95)

#Table of means and sd of political interest per field of study/work

group_by(data_pooled, Politics) %>%
  summarise(
    count = n(),
    mean = mean(polinterest, na.rm = TRUE),
    sd = sd(polinterest, na.rm = TRUE)
  )

#Table of means and sd of political interest per gender
group_by(data_pooled, gender)%>%
  summarise(
    count = n(),
    mean = mean(polinterest, na.rm = TRUE),
    sd = sd(polinterest, na.rm = TRUE)
  )

# Correlation between Political interest and political knowledge
cor.test(data_pooled$polinterest, data_pooled$polknowledge)

#
apa.aov.table(
aov,
"sample2_ex_interest.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)
```

### Recoding the Politics group

Table S16
```{r Field of study/work not including History stage 2}
#In the original study, the Politics group was composed by bachelor students of Political Science, Public Administration, Sociology and Governance. 
data_pooled$Politics_new <- ifelse(grepl("1", data_pooled$Field_1) == TRUE  |
                          #grepl("1", data_pooled$Field_5) == TRUE  |
                          grepl("1", data_pooled$Field_11) == TRUE  |
                          grepl("1", data_pooled$Field_12) == TRUE  |
                          grepl("1", data_pooled$Field_27) == TRUE, "Yes", "No")

ancova <- aov(polknowledge ~ polinterest + gender*Politics_new*condition, data_pooled)
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

sjstats::anova_stats(ancovafit)
sjstats::eta_sq(ancovafit, partial = T, ci.lvl = .95)

apa.aov.table(
ancova,
"sample2_ex_nohist.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)
```

### Removing field of study/work from analysis

Table S17
```{r Not including field of study/work stage 2}
ancova <- aov(polknowledge ~ polinterest + gender*condition, data_pooled)
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

sjstats::anova_stats(ancovafit)
sjstats::eta_sq(ancovafit, partial = T, ci.lvl = .95)

apa.aov.table(
ancova,
"sample2_ex_nofield.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)
```

### Participants’ ranking of Politics as an important topic for their studies/work

Table S18
```{r Topic ranking as covariate stage 2}
# relabel topics ranking with topic name
names(data_pooled)[217:243] <- c("rank_politics", 
                                     "rank_psychology",
                                     "rank_economics",
                                     "rank_biology",
                                     "rank_history",
                                     "rank_geography",
                                     "rank_chemistry",
                                     "rank_physics",
                                     "rank_business",
                                     "rank_foodindustry",
                                     "rank_sociology",
                                     "rank_publicadministration",
                                     "rank_mathematics",
                                     "rank_law",
                                     "rank_accounting",
                                     "rank_engineering",
                                     "rank_pedagogy",
                                     "rank_education",
                                     "rank_healthcare",
                                     "rank_socialservices",
                                     "rank_arts",
                                     "rank_finance",
                                     "rank_transportation",
                                     "rank_retail",
                                     "rank_it",
                                     "rank_realestate",
                                     "rank_government") 
#
#table(data_pooled$Politics, data_pooled$rank_politics, useNA = "always")
# reverse ranking so that NA = 0, a rank of 1 = 4 --> ordinal scale
data_pooled$politics_ranking <- ifelse(is.na(data_pooled$rank_politics), 0,
                                       ifelse(data_pooled$rank_politics == 1, 4,
                                              ifelse(data_pooled$rank_politics == 2, 3,
                                                     ifelse(data_pooled$rank_politics == 3, 2, 1))))
#
#table(data_pooled$politics_ranking, useNA = "always")
#
ancova <- aov(polknowledge ~ polinterest + politics_ranking + gender*Politics*condition, data_pooled)
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

sjstats::eta_sq(ancovafit, partial = T, ci.lvl = .95)

apa.aov.table(
ancova,
"sample2_ex_rank.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)
```

### Participants’ attributed importance to Politics

Table S19
```{r Topic importance as covariate and between-subjects factor stage 2}
# relabel topics importance with topic name
names(data_pooled)[249:275] <- c("importance_politics", 
                                     "importance_psychology",
                                     "importance_economics",
                                     "importance_biology",
                                     "importance_history",
                                     "importance_geography",
                                     "importance_chemistry",
                                     "importance_physics",
                                     "importance_business",
                                     "importance_foodindustry",
                                     "importance_sociology",
                                     "importance_publicadministration",
                                     "importance_mathematics",
                                     "importance_law",
                                     "importance_accounting",
                                     "importance_engineering",
                                     "importance_pedagogy",
                                     "importance_education",
                                     "importance_healthcare",
                                     "importance_socialservices",
                                     "importance_arts",
                                     "importance_finance",
                                     "importance_transportation",
                                     "importance_retail",
                                     "importance_it",
                                     "importance_realestate",
                                     "importance_government") 
#
#
# as covariate
#
ancova <- aov(polknowledge ~ polinterest + importance_politics + gender*Politics*condition, data_pooled)
ancovafit <- car::Anova(ancova, type = 3)
ancovafit
#
sjstats::eta_sq(ancovafit, partial = T, ci.lvl = .95)
#
apa.aov.table(
ancova,
"sample2_ex_importance.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)
```

### Sample composition

Table S20
```{r younger sample stage 2}
# Ihme& Tausendpfund (2018, p. 44)
# Test subjects included N = 377 students of psychology (n = 236) and politics (n = 141) 
# at the FernUniversität in Hagen (M = 31 years, SD = 12.04; 58.9% female).

# Our sample:
# 1502 participants (Mage = 45.87 years, SDage = 17.35, 48.74% female)

psych::describe(data_pooled$Age)
#
data_pooled$age_group <- ifelse(data_pooled$Age < 44, "young", 
                                ifelse(data_pooled$Age > 43, "old", NA))
table(data_pooled$age_group, useNA = "ifany")
#
psych::describeBy(data_pooled$Age, group = data_pooled$age_group) # M = 31, SD = 7.17
#
table(data_pooled$age_group, data_pooled$Education)
#
data_pooled$age_education <- ifelse(data_pooled$age_group == "young" & data_pooled$Education >= 3, "young and higher ed",
                                    ifelse(data_pooled$age_group == "young" & data_pooled$Education < 3, "yound and lower ed",
                                           ifelse(data_pooled$age_group == "old" & data_pooled$Education >= 3, "old and higher ed",
                                                  ifelse(data_pooled$age_group == "old" & data_pooled$Education < 3, "old and lower ed", NA))))
#
table(data_pooled$age_education, useNA = "ifany") # 522 participants young and higher ed
#
data_younged <- data_pooled[which(data_pooled$age_education == "young and higher ed"),]

ancova <- aov(polknowledge ~ polinterest + gender*Politics*condition, data_younged)
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

#sjstats::anova_stats(ancovafit)
sjstats::eta_sq(ancovafit, partial = T, ci.lvl = .95)

apa.aov.table(
ancova,
"rr_newsample.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)
#
#ancova <- aov(polknowledge ~ polinterest + gender*Politics*condition*age_education, data_pooled)
#emmeans.ancova <- emmeans(ancova, specs = c("age_education"))
#emmeans.ancova 

#contrast <- contrast(emmeans.ancova, "pairwise", 
 #                       simple = "gender", # "each" gives all contrasts
  #                      combine = T, 
   #                     adjust = "bonferroni")

#contrast
```

### Measurement Invariance

Table S21
```{r Measurement Invariance}
PolKnow_items_all_MI <- PolKnow_items_all
names(PolKnow_items_all_MI) <- c("i1", "i2", "i3", "i4", "i5", "i6", "i7", "i8", "i9", "i10", "i11", "i12", "i13", "i14", "i15", "i16", "i17", "i18", "i19", "i20")
PolKnow_items_all_MI$gender  <- data_pooled$gender
PolKnow_items_all_MI$genderN <- as.numeric(PolKnow_items_all_MI$gender)-1

# Measurement Invariance
# ----------------------

library(difR)
library(ltm)

Lord.dif.1PL <- difLord(PolKnow_items_all_MI[, 1:20], group = PolKnow_items_all_MI[, 22], focal.name = 1, model = "1PL")
plot(Lord.dif.1PL)

Lord.dif.2PL <- difLord(PolKnow_items_all_MI[, 1:20], group = PolKnow_items_all_MI[, 22], focal.name = 1, model = "2PL")
plot(Lord.dif.2PL)

# 3PL doesn't estimate
#Lord.dif.3PL <- difLord(PolKnow_items_all_MI[, 1:20], group = PolKnow_items_all_MI[, 22], focal.name = 1, model = "3PL", discr = NULL)

library("mirt")

#MIRT
mod_configural <- multipleGroup(PolKnow_items_all_MI[, 1:20], 1, group = PolKnow_items_all_MI[, 21], SE = TRUE)

M2(mod_configural) # OK ish
plot(mod_configural)

mirtCluster()
DRF(mod_configural, draws=500, plot=TRUE)

# constrain slopes within each group to be equal (but not across groups)
model2 <- 'F1 = 1-20
           CONSTRAIN = (1-20, a1)'
mod_configural2 <- multipleGroup(PolKnow_items_all_MI[, 1:20], model2, group = PolKnow_items_all_MI[, 21], SE = TRUE)

summary(mod_configural2)
DRF(mod_configural2, draws=500, plot=TRUE)


DTF(mod_configural, draws = 1000) 
DTF(mod_configural, draws = 1000, plot='sDTF')
```

## Analyses of variance with mean-centered covariates and orthogonal contrasts

Table S22
```{r orthogonal factors ANOVA stage 1b}
# Given that type III analyses only give sensible and informative results when covariates are mean-centered and factors are coded with orthogonal contrasts, we adjusted both, covariates and factor variables.

aov <- aov(polknowledge ~ gender*Politics*condition, data_stage1, contrasts=list(gender=contr.sum, Politics=contr.sum, condition=contr.sum))
aovfit <- car::Anova(aov, digits = 3, type = 3)
aovfit

sjstats::eta_sq(aovfit, partial = T, ci.lvl = .95)

apa.aov.table(
aov,
"sample1_h1_corrected.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)

```

Table S23
```{r mean-centered covariate and orthogonal factors ANCOVA stage 1}
# Given that type III analyses only give sensible and informative results when covariates are mean-centered and factors are coded with orthogonal contrasts, we adjusted both, covariates and factor variables.
data_stage1$polinterest_c <- center_scale(data_stage1$polinterest)

ancova <- aov(polknowledge ~ polinterest_c + gender*Politics*condition, data_stage1, contrasts=list(gender=contr.sum, Politics=contr.sum, condition=contr.sum))
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

#sjstats::anova_stats(ancovafit)
sjstats::eta_sq(ancovafit, partial = T, ci.lvl = .95)

apa.aov.table(
ancova,
"sample1_h2_corrected.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)
```

Table S24
```{r orthogonal factors ANOVA stage 2}
# Given that type III analyses only give sensible and informative results when covariates are mean-centered and factors are coded with orthogonal contrasts, we adjusted both, covariates and factor variables.

aov <- aov(polknowledge ~ gender*Politics*condition, data_pooled, contrasts=list(gender=contr.sum, Politics=contr.sum, condition=contr.sum))
aovfit <- car::Anova(aov, digits = 3, type = 3)
aovfit

#sjstats::anova_stats(aovfit, digits = 3, partial = T)
sjstats::eta_sq(aovfit, partial = T, ci.lvl = .95)

apa.aov.table(
aov,
"sample2_h1_corrected.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)

```

Table S25
```{r mean-centered covariate and orthogonal factors ANCOVA stage 2b}
# Given that type III analyses only give sensible and informative results when covariates are mean-centered and factors are coded with orthogonal contrasts, we adjusted both, covariates and factor variables.
data_pooled$polinterest_c <- center_scale(data_pooled$polinterest)

ancova <- aov(polknowledge ~ polinterest_c + gender*Politics*condition, data_pooled, contrasts=list(gender=contr.sum, Politics=contr.sum, condition=contr.sum))
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

#sjstats::anova_stats(ancovafit)
sjstats::eta_sq(ancovafit, partial = T, ci.lvl = .95)

apa.aov.table(
ancova,
"sample2_h2_corrected.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)

```

Table S26
```{r Dont know answers mean-centered covariate and orthogonal factors ANCOVA stage 2}
# Given that type III analyses only give sensible and informative results when covariates are mean-centered and factors are coded with orthogonal contrasts, we adjusted both, covariates and factor variables.

ancova <- aov(PKT_DK ~ polinterest_c + gender*condition, data_pooled, contrasts=list(gender=contr.sum, condition=contr.sum))
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

apa.aov.table(
ancova,
"sample2_ex_DK_corrected.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)


```

Table S27
```{r Incorrect answers mean-centered covariate and orthogonal factors ANCOVA stage 2}
# Given that type III analyses only give sensible and informative results when covariates are mean-centered and factors are coded with orthogonal contrasts, we adjusted both, covariates and factor variables.
ancova <- aov(PKT_incorrect ~ polinterest_c + gender*condition, data_pooled, contrasts=list(gender=contr.sum, condition=contr.sum))
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

apa.aov.table(
ancova,
"sample2_ex_inc_corrected.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)

```

Table S28
```{r DK/incorrect proportion mean-centered covariate and orthogonal factors ANCOVA stage 2}
# Given that type III analyses only give sensible and informative results when covariates are mean-centered and factors are coded with orthogonal contrasts, we adjusted both, covariates and factor variables.

#data_pooled$polinterest_c <- center_scale(data_pooled$polinterest)
data_pooled$prop.dk_c <- center_scale(data_pooled$prop.dk)

ancova <- aov(polknowledge ~ polinterest_c + prop.dk_c + gender*Politics*condition, data_pooled, contrasts=list(gender=contr.sum, Politics=contr.sum, condition=contr.sum))
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

apa.aov.table(
ancova,
"sample2_ex_prop_corrected.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)
```

Table S29
```{r Attempted items mean-centered covariate and orthogonal factors ANCOVA stage 2}
# Given that type III analyses only give sensible and informative results when covariates are mean-centered and factors are coded with orthogonal contrasts, we adjusted both, covariates and factor variables.

ancova <- aov(attempt ~ polinterest_c + gender*condition, data_pooled, contrasts=list(gender=contr.sum, condition=contr.sum))
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

apa.aov.table(
ancova,
"sample2_ex_attempt_corrected.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)
```

Table S30
```{r Political Interest orthogonal factors ANOVA stage 2}
# Given that type III analyses only give sensible and informative results when covariates are mean-centered and factors are coded with orthogonal contrasts, we adjusted both, covariates and factor variables.
aov <- aov(polinterest ~ gender*Politics*condition, data_pooled, contrasts=list(gender=contr.sum, Politics=contr.sum, condition=contr.sum))
aovfit <- car::Anova(aov, type = 3)
aovfit

apa.aov.table(
aov,
"sample2_ex_interest_corrected.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)
```

Table S31
```{r Field of study/work not including History mean-centered covariate and orthogonal factors ANCOVA stage 2}
# Given that type III analyses only give sensible and informative results when covariates are mean-centered and factors are coded with orthogonal contrasts, we adjusted both, covariates and factor variables.

ancova <- aov(polknowledge ~ polinterest_c + gender*Politics_new*condition, data_pooled, contrasts=list(gender=contr.sum,Politics_new=contr.sum, condition=contr.sum))
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

apa.aov.table(
ancova,
"sample2_ex_nohist_corrected.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)

```

Table S32
```{r Not including field of study/work mean-centered covariate and orthogonal factors ANCOVA stage 2}
# Given that type III analyses only give sensible and informative results when covariates are mean-centered and factors are coded with orthogonal contrasts, we adjusted both, covariates and factor variables.
ancova <- aov(polknowledge ~ polinterest_c + gender*condition, data_pooled, contrasts=list(gender=contr.sum, condition=contr.sum))
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

apa.aov.table(
ancova,
"sample2_ex_nofield_corrected.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)

```

Table S33
```{r Topic ranking as covariate mean-centered covariate and orthogonal factors ANCOVA stage 2}
# Given that type III analyses only give sensible and informative results when covariates are mean-centered and factors are coded with orthogonal contrasts, we adjusted both, covariates and factor variables.

data_pooled$politics_ranking_c <- center_scale(data_pooled$politics_ranking)

ancova <- aov(polknowledge ~ polinterest_c + politics_ranking_c + gender*Politics*condition, data_pooled, contrasts=list(gender=contr.sum, Politics=contr.sum, condition=contr.sum))
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

apa.aov.table(
ancova,
"sample2_ex_rank_corrected.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)

```

Table S34
```{r Topic importance as covariate and between-subjects factor mean-centered covariate and orthogonal factors ANCOVA stage 2}
# Given that type III analyses only give sensible and informative results when covariates are mean-centered and factors are coded with orthogonal contrasts, we adjusted both, covariates and factor variables.

data_pooled$importance_politics_c <- center_scale(data_pooled$importance_politics)

ancova <- aov(polknowledge ~ polinterest_c + importance_politics_c + gender*Politics*condition, data_pooled, contrasts=list(gender=contr.sum, Politics=contr.sum, condition=contr.sum))
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

apa.aov.table(
ancova,
"sample2_ex_importance_corrected.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)
```

Table S35
```{r younger sample mean-centered covariate and orthogonal factors ANCOVA stage 2}
# Given that type III analyses only give sensible and informative results when covariates are mean-centered and factors are coded with orthogonal contrasts, we adjusted both, covariates and factor variables.

ancova <- aov(polknowledge ~ polinterest_c + gender*Politics*condition, data_younged, contrasts=list(gender=contr.sum, Politics=contr.sum, condition=contr.sum))
ancovafit <- car::Anova(ancova, type = 3)
ancovafit

apa.aov.table(
ancova,
"sample2_rr_newsample_corrected.doc",
table.number = NA,
conf.level = 0.95,
type = 3
)
```

## Reporting Standards

```{r average completion time stage 1 and stage 2}
psych::describe(data_stage1$Duration__in_seconds_)
psych::describe(data_pooled$Duration__in_seconds_)
```

## Saving all data 

We are saving data in .RData format that is produced at the end of the knitted report so that replicators have access to all data sources and outputs generated by our analyses. With this data one should be able to know what results we found with our analyses (thereby hopefully circumventing computational reproducibility issues).

```{r}
save.image("Reproducible Data for all analyses.RData")
```
